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Biochemical networks are the analog computers of life. They allow living cells to control a large 
number of biological processes, such as gene expression and cell signalling. In biochemical networks, 
the concentrations of the components are often low. This means that the discrete nature of the 
reactants and the stochastic character of their interactions have to be taken into account. Moreover, 
the spatial distribution of the components can be of crucial importance. However, the current 
numerical techniques for simulating biochemical networks either ignore the particulate nature of 
matter or treat the spatial fluctuations in a mean-field manner. We have developed a new technique, 
called Green's Function Reaction Dynamics (GFRD), that makes it possible to simulate biochemical 
networks at the particle level and in both time and space. In this scheme, a maximum time step is 
chosen such that only single particles or pairs of particles have to be considered. For these particles, 
the Smoluchowski equation can be solved analytically using Green's functions. The main idea of 
GFRD is to exploit the exact solution of the Smoluchoswki equation to set up an event-driven 
algorithm. This allows GFRD to make large jumps in time when the particles are far apart from 
each other. Here, we apply the technique to a simple model of gene expression. The simulations 
reveal that the scheme is highly efficient. Under biologically relevant conditions, GFRD is up to six 
orders of magnitude faster than conventional particle-based techniques for simulating biochemical 
networks in time and space. 

PACS numbers: 87.16.Yc,87.16.Ac,2.70.-c,5.40.-a 



I. INTRODUCTION 

Organisms can be viewed as information processing 
machines. Even relatively simple organisms, such as the 
bacterium Escherichia coli, can perform fairly complex 
computations such as: if lactose is present and 
not glucose is present, then use lactose. Re- 
cent technological developments have made it possible 
to obtain information on the regulatory architecture of 
the cell on an unprecedented scale. In addition, exten- 
sive databases are now available that catalog biochemical 
pathways. Nevertheless, our understanding of the mech- 
anisms that allow living cells to process information is 
still limited. One important reason for this is that these 
mechanisms are controlled by stochastic processes. 

In the living cell, computations are performed by 
molecules that chemically and physically interact with 
each other. These components, that form what is called 
a biochemical network, behave stochastically. They often 
move in an erratic fashion, namely by diffusion, and act 
upon each other in a stochastic manner - chemical re- 
actions, and equally important, physical interactions are 
probabilistic in nature. These factors become particu- 
larly important when the concentrations are low. In the 
living cell, this is often the case, and, as a result, bio- 
chemical networks can be highly stochastic 0, Q . In this 
respect, it is a remarkable fact that many biological pro- 
cesses operate reliably with surprisingly small numbers 
of molecules. 



Another important reason for our limited understand- 
ing of biochemical networks is that they operate not only 
in time, but also in space. In the living cell, signals of- 
ten have to be transmitted from one place to the next by 
the diffusive motion of "messenger" molecules; their con- 
centrations can be non-uniform, and more importantly, 
their low mobility can limit the response time of the net- 
work. Moreover, many processes, such as the immune 
response, involve a complex spatial reorganization of the 
reactants. Finally, a large number of biological activities 
require the localized assembly of a complex of proteins. 
All these processes can only be accurately investigated 
using techniques that resolve the network in time and 
space. 

In principle, computer simulations are ideally suited 
for studying how biochemical networks reliably process 
information in time and space. However, the current nu- 
merical techniques are of limited use for this purpose. Ta- 
bic [I] gives an overview of the commonly used techniques 
for analysing biochemical networks. The conventional 
approach is to write down the macroscopic rate equations 
and to solve the corresponding differential equations nu- 
merically. In this method, the evolution of the network is 
deterministic. It is implicitly assumed that the concen- 
trations are large and that fluctuations can be neglected. 
The effect of fluctuations is often included by adding a 
noise term to the macroscopic rate equations 3J. How- 
ever, at low concentrations, this approach is bound to 
fail, as demonstrated by Togashi and Kaneko and 



2 



Shnerb and coworkers |5(. At low concentrations, we 
have to recognize the discrete nature of the reactants 
and the stochastic character of their interactions. Cur- 
rently, two techniques exist for simulating biochemical 
networks at the particle level |g, [7| • Both of these are 
consistent with the chemical master equation. However, 
the chemical master equation relies on the assumption 
that there are many non-reactive collisions to stir the 
system between the reactive collisions. In effect, it is im- 
plicitly assumed that at each instant the particles are uni- 
formly distributed in space. This is a serious limitation. 
As discussed above, the functioning of many biochemical 
networks depends crucially on their spatial organization. 
Moreover, fluctuations of the components in space can 
be a major source of noise in biochemical networks. 

Clearly, in many cases we have to describe biochemical 
networks both temporally and spatially at the molecu- 
lar level. However, the current techniques for simulating 
biochemical networks cannot accomplish this since they 
either ignore the particulate nature of matter, or treat 
the spatial fluctuations in a mean-field manner. We have 
developed a new technique, named Green's Function Re- 
action Dynamics (GFRD), that makes it possible to sim- 
ulate biochemical networks at the particle level and in 
time and space. This technique, which we describe in 
detail in the next section, is an event-driven algorithm 
in which the particles are propagated in space between 
reaction events using Green's functions. In the subse- 
quent section, we apply GFRD to a simple model of gene 
expression. The calculations show that the technique is 
highly efficient. We thus believe that GFRD brings sim- 
ulating biochemical networks at the particle level and in 
time and space within reach. 



II. NUMERICAL TECHNIQUE 
A. Introduction 

In this section we discuss the main ideas underlying 
Green's Function Reaction Dynamics (GFRD). GFRD is 
a stochastic scheme that combines the propagation of the 
particles in space with the reactions between them in an 
efficient manner. Most proteins are believed to be trans- 
ported by their diffusive motion, although in eukaryotic 
cells, cytoskeletal networks and motor proteins may fa- 
cilitate the transport of molecules 8]. Such mechanisms 
of active transport have not been observed in bacteria, 
where diffusion is believed to be the primary means of 
intracellular movement. For now, we will assume that 
the network components are transported by diffusion, al- 
though the technique could be extended to include active 
transport as well. 

Two approaches seem to be potentially useful for sim- 
ulating biochemical networks at the particle level and in 
time and space. The first is to let the particles undergo a 
random walk on a lattice and to let reaction partners re- 
act with a certain probability when they happen to meet 



each other. This technique has a number of limitations, 
the most important of which are that the physical dimen- 
sions of the particles and the interactions between them 
cannot conveniently be described. 

Brownian Dynamics is a more appealing technique. 
This is a stochastic dynamics scheme, in which the parti- 
cles are propagated in space according to the overdamped 
limit of the Langevin equation. In Brownian Dynamics, 
the solvent is considered implicitly; only the solute par- 
ticles are considered explicitly. The forces experienced 
by these particles contain two parts: a conservative part, 
which arises from the interactions with the other solute 
particles, and a random part. The latter is the dynamical 
remnant of the solvent - the solutes are thought to experi- 
ence random forces from the solvent. Via the fluctuation- 
dissipation theorem and the Einstein relation, the ran- 
dom forces are related to the diffusion constant of the 
particles. To be more explicit, the equations of motion 
for the solute particles are given by: 

*s = ^(F s + 5F s ). (1) 
k B l 

Here, r s denotes the position of solute particle s, D s is 
the diffusion constant of solute particle s, k B T is Boltz- 
mann's constant times temperature, F s is the conserva- 
tive force exerted by the other solute particles, and SF S 
is the random force that arises from the interactions with 
the solvent. 

Brownian Dynamics has a number of advantages over 
lattice-based techniques: the particles move in contin- 
uum space; the interactions between particles - the po- 
tential of mean force - can easily be described; excluded 
volume effects are taken into account naturally; and a 
different diffusion constant can be assigned to each type 
of particles. 

In principle, chemical reactions can be straight- 
forwardly implemented into the Brownian Dynamics 
scheme: the particles are propagated according to Eq. 
and when two reaction partners happen to meet each 
other, they can react with a probability that is consistent 
with the rate constant. However, the major drawback of 
such a scheme is that very small time steps are needed in 
order to resolve the collision events. This makes brute- 
force Brownian Dynamics a very inefficient scheme to 
simulate biochemical networks. 

The main idea of GFRD is to combine in one step 
the propagation of the particles in space with the re- 
actions between them. To see how this can be accom- 
plished, one should realize that Brownian Dynamics is, 
in effect, a numerical procedure for solving the Smolu- 
chowski equation 0. However, for a pair of particles, 
that not only move diffusively, but also can react ac- 
cording to A + B — > C + D + the Smoluchowski 
equation can be solved analytically using Green's func- 
tions. The Green's function for the pair of particles A 
and B, p(r, i|r , to), yields the probability that the inter- 
particle vector ro at time to becomes r at a later time t. 
The essence of GFRD is to exploit this exact solution for 
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Description 




Accounts for spatial 
extent of network 


Incorporates 
fluctuations 




Ordinary differential 
equations 


No 


No 


Continuum 


Stochastic 

differential equations 


No 


Only at high 
concentrations 




Reaction diffusion 
equations 


Yes 


1NO 


Particle-based 


Chemical master 
equation 


No 


Yes 


GFRD 


Yes 


Yes 



TABLE I: Overview of the commonly used techniques and the newly developed technique, called Green's Function Reaction 
Dynamics (GFRD), to simulate biochemical networks. GFRD takes both the discrete nature and the spatial distribution of 
the reactants into account. 



a pair of particles to set up an event-driven algorithm. 
This allows GFRD to make large jumps in time when the 
particles are far apart from each other. In biochemical 
networks this is often the case as the reactant concentra- 
tions are usually low. GFRD is therefore ideally suited 
for biochemical networks. Below we describe the scheme 
in more detail. 



C. Monomolecular reactions - the Green's function 
for single particles 

In this section, we consider the propagation of a single 
particle. We assume that the particle is spherical in shape 
and moves by free diffusion with a diffusion constant D. 
The diffusive motion of the particle is described by the 
Einstein diffusion equation 

d tPl {r,t\r Q ,t ) = DV 2 Pl {r,t\r ,t ). (2) 



B. Overview of the algorithm 



Let us imagine a configuration of reactants as shown 
in Fig. n The circles indicate the distances the parti- 
cles can travel in a time step At. For a particle that 
moves by free diffusion with a diffusion constant D, that 
distance scales as (|Ar| 2 ) ~ DAt. Intersecting circles 
represent particles that may meet within the time step 
At. Isolated particles and pairs of interacting particles 
can be propagated analytically using Green's functions, 
as we will discuss in detail in the next section. Clearly, 
the larger the time step, the larger the circles and the 
greater the probability that reaction partners meet and 
react with each other. However, we cannot make the time 
step arbitrarily large: if a given particle can collide with 
more than one other particle during a time step, then 
propagating the particles becomes a many-body problem 
that we can not solve analytically. The size of a time 
step is thus limited by the requirement that each parti- 
cle can interact with at most one other particle during 
that time step. This constraint sets an upper limit on 
the magnitude of a time step in our algorithm; we will 
call it Ai max . However, provided that we consider times 
At smaller than Ai max , the problem can be reduced to 
that of propagating single particles and pairs of particles. 
This can be solved analytically using Green's functions, 
as we will now describe. 




FIG. 1: Determination of the maximum time step, A£ max . 
The maximum size of the time step is set by the requirement 
that each particle can interact with at most one other par- 
ticle during that time step; each particle i can thus travel 
a distance of at most Ar maXj i during a time step, as indi- 
cated by the circles. We have used the operational criterion 
Ar maXj i = H-y/6DiAt ma x,i, with Di being the diffusion con- 
stant of particle i. A value of H = 3 was found to yield a 
good conservation of the correct steady-state distribution. In 
this example, each particle is assumed to have the same diffu- 
sion constant; the time step is limited by the constraint that 
particles i and k should not interact as particle i can already 
interact with particle j. Note that with this maximum time 
step the many-body problem of propagating the N particles 
is reduced to that of propagating single particles and pairs of 
particles. 
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Here, pi(r,t|ro,to) is the probability that the particle is 
at position r at time t, given that is was at r at time 
to. This diffusion equation can be solved subject to the 
initial condition and the boundary condition 

Pi(r,t |r ,to) = S(t-t ), (3) 
pi(|r| -+oo,i|r ,i ) = 0, (4) 

respectively. The solution pi(r,t\r ,t ) is known as a 
Green's function. It is given by the well-known expres- 
sion 



pi(r,i|r ,to) 



1 



[47TD(t - t )f 2 



exp 



4D(t-to) 



( 5 ) 

We now consider the case in which the particle does not 
only move diffusively, but also can "decay" according to 



A^B + C. 



(6) 



We will assume that, if the reaction happens, it happens 
instantaneously. This means that the reaction can be de- 
coupled from the diffusive motion of the particle. If the 
reaction is a Poisson process with k d dt being the proba- 
bility that a reaction occurs in an infinitesimal time inter- 
val dt, then the probability that the next reaction occurs 
between t and t + dt, is 



q d (t\t )dt = k d exp [-k d (t - t )} dt. 



(7) 



In section III El we will use Eqs. |S] and [7| to set up an 

event-driven algorithm. 



D. Bimolecular reactions - the Green's function for 
pairs of particles 

In this section, we consider one pair of particles A and 
B that can react according to 



A + B 



C + D.. 



(8) 



We again assume that the particles A and B are spherical 
in shape and move by their diffusive motion; the diffu- 
sion constants for particle A and B are Da and Db, 
respectively. Furthermore, we assume that the parti- 
cles react with an intrinsic rate constant k a when they 
have approached each other within the reaction distance 
a = [d,A + ds)/2, where dA and ds are the diameters of 
particles A and B, respectively. In addition, we imagine 
that the particles interact with each other via a potential 
U(r), where r = yb — ya- The force acting on particle B 
is thus given by — V bU(y) — F(r), while the force acting 
on particle A is given by — F(r). 

We aim to derive the distribution function 
P2 (ta , rfl , 1 1 r , rso , to ) , which yields the probabil- 
ity that the particles A and B are at positions ya and 
yb at time t, given that they were at r^o and r^o at 



time to. This distribution function satisfies for |r| > a 
the following Smoluchowski equation [|| 

dtP2(rA,Y B ,t\Y A 0,rBO,to) = 

[D A V 2 A + D B V 2 B - DbPVb ■ F(r) + D A f3V A • F(r)] 
x p 2 {y a , r B , t\Y A o, r B o, to). (9) 

It will be convenient to make a coordinate transformation 
R = ^Db/Daya + ^Da/DbYb, 



Y = Y B - Y A , 

and to define the operators 

V R - d/dR, 
V r = d/dr. 

Eq. E|can then be rewritten as: 

d t p 2 (R,Y,t\R ,ro,t ) = ( J D j4 + J Ds)[V^ + V r -(V r 

X p 2 (R, Y,t\R ,Y ,t Q ), 

Irl > a. 



(10) 
(11) 

(12) 
(13) 

-F(r))] 
(14) 



It is seen that Eq. ^] describes two independent random 
processes - free diffusion in the coordinate R and diffu- 
sion with a drift in the coordinate r. This means that the 
distribution function p 2 {ya, ?b, t\YAo, ^BOi to) can be fac- 
torized asp^(R, t\Ro, to)p\{Y, t|r , to) and that the above 
equation can be reduced to one Smoluchowski equation 
for the coordinate R and one for the coordinate r: 

%^(R,t|R ,io) = (D A + D B )V^ 

xpf(R,t\R ,to), (15) 
d t p T 2 (Y,t\Y ,t ) = {D A +D B )V r -(V r -F(r)) 

xp r 2 (r,t\r,t ),\r\><7- (16) 

Eqn. El describes the free diffusive motion of the coor- 
dinate R. The solution of that equation, for the initial 
condition p 2 '(R, to\Ro, to) = S(R — Ro) and boundary 
condition p^(|R| — > oo, t|Ro,to) = 0, is 



pf(R,t\R ,t ) 



exp 



[R-RoT 



A{D A +D B ){t-t ) 



[4TT(DA + D B )(t-to)} 



3/2 ' 



(17) 



The non-trivial solution is that of the Smoluchowski 
equation for the inter-particle vector r. This solution 
also has to take into account the reactions between A 
and B. We will incorporate the reaction as a boundary 
condition on the solution of the Smoluchoswki equation. 
To be more explicit, the initial condition and boundary 
conditions for the coordinate r are given by 



P2( r >*o|ro,*o) = 6(t- to), 
pl(\r\ -> oo,t\r ,t ) = 0, 

-j(a,t\r ,t ) = 4:TTa 2 D 



(18) 
(19) 

d_ 

Or 

x P r 2 (r,t\Y ,t )\ lr \ =<T , 
k a p r 2 (\r\ =o-,t|r ,to), (20) 



F(r) 
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where d/d r denotes a derivative with respect to the inter- 
particle separation r. It is seen that the reaction enters 
the problem as a third boundary condition on the so- 
lution of the Smoluchowski equation. Here j(er, t|ro, to) 
is the outward radial flux of probability p^i^, t\vo, to) 
through the "contact" surface of area 47r<r 2 . The bound- 
ary condition, also known as a radiation boundary con- 
dition [llj, [l2(, states that this radial flux of probabil- 
ity equals the intrinsic rate constant times the proba- 
bility that the particles A and B are, in fact, in con- 
tact. In the limit k a — > oo, the radiation boundary 
condition reduces to an absorbing boundary condition 
j?2 ( I r I — 0"> ^l r 0: £o) = 0, while in the limit k a — ► the ra- 
diation boundary condition reduces to a reflecting bound- 
ary condition. 

The Green's function p\(r, t|ro, to) is derived in ap- 
pendix A for the case in which F(r) = for |r| > a; for 
cases in which F(r) ^ for |r| > a, the Green's func- 
tions could, depending upon the interaction potential, 
either be obtained analytically or numerically. Here, we 
will discuss some useful quantities that can be derived 
from the Green's function. The first quantity of interest 
is the probability that a pair of particles, initially sepa- 
rated by ro, survives and does not recombine by time t. 
This so-called survival probability is given by 

S a (t\r ,t ) = / drp5(r,t|r ,to). (21) 

^ |rj >cr 

Clearly, 5 a (0|r ,t ) = 1 for |r | > a. The second quan- 
tity of interest is the reaction rate q a (t\ro, to), which is 
defined as the probability per unit time that a pair, ini- 
tially separated by ro, reacts at time t. It is related to 
the survival probability by 

u\ + \ dS a (t\r ,to) 
q a {t\r , t ) = — . (22) 

Since the reactions are assumed to occur only at contact, 
the reaction rate is also given by the flux at contact 

q a (t\r ,t ) = -j(a,t\r ,t ). (23) 

The above equation can also be obtained from Ea. 1 161 and 
Eq. 1211 and by using the fact that the flux at |r| — > oo 
vanishes. 

The reaction rate g a (t|i"o,io) can be interpreted as the 
probability that the next reaction for a pair of particles, 
initially separated by ro, occurs at time t. This is used to 
set up the GFRD event-driven algorithm, which we will 
describe in the next section. 



E. Outline of the algorithm 

To explain the essence of the algorithm, it will be in- 
structive to consider a single particle of type A that can 
react with a single particle of type B according to the 
following scheme 

A + B^C. (24) 



Furthermore, it will be useful to imagine that these par- 
ticles are surrounded by neighbors, the presence of which 
limit the size of the time step to At max . As a function 
of time, the system will flip-flop between the associated 
state C and the dissociated state A + B. The GFRD 
event-driven algorithm to propagate this system would 
consist of iterating the following steps: 

1. If the system is in the dissociated state A + B 7 
then draw a next association time t according to 
q a {t\ro,t ) (Eq.E). 

a) If (t — to) > At max , then the two particles 
will not react within the time step; new positions 
for A and B at time to + At max are obtained from 
pf{R,t + At roax |R ,t ) (Eq. H3 and p r 2 (r,t + 
At max |r ,to) (Ea. I37|l with R and r as given by 
Eas.HlHandlTTl 

b) If [t — to) < At max , then the next reaction 
will occur within the time step; a new position for 
particle C at time t is obtained fromp^-(R, t|R , to) 
(Eq.HZJ. 

2. If the system is in the associated state C, then draw 
a next dissociation time from qd(t\to) (Eq. 0). 

a) If (t — t ) > At max , then particle C will 
not have decayed by to + At max ; a new position for 
particle C, re, at time to 4- At max is obtained from 
Pi(rc»*o + At max |rco,io) (Eq.EJ); 

b) If (t — to) < At max , the next reaction will 
occur within the maximum time step; the parti- 
cles A and B are placed at time t adjacent to 
each other at positions around vq as obtained from 
Pi(rc,t|rco,*o) (see Eq.|S}. 

The procedure outlined above forms the heart of the 
algorithm. The crux of the method is to choose the 
maximum time step such that only monomolecular re- 
actions or bimolecular reactions have to be considered. 
This makes it possible to use the exact solution of the 
Smoluchowski equation to propagate the system to the 
next reaction event in a single step. The full algorithm 
for a system of N particles thus consists of iterating the 
following steps: 

1. Determine maximum time step At max . The maxi- 
mum time step is determined by the condition that 
only single particles or pairs of particles have to be 
considered (see section III Bl and Fig. QJ. For each 
particle i, we determine the maximum time step 
At max .i, such that it can interact with at most one 
other particle. The maximum global time step is 
then given by 

At max = min({At maXi J). (25) 

In order to determine At max i for particle i, we as- 
sume that during that step the particle can travel 
at most a distance Ar max .i = HyJ6DiAt, where Di 
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is the diffusion constant of particle i. We find that 
H = 3 suffices to preserve the correct steady-state 
distribution. 

2. Determine next reaction and next reaction time . 
We first construct a list of possible reactions {R u }- 
With each reaction R v , we associate a survival 
probability function S v (t — to) and a next-reaction 
distribution function q v (t — to)', the two are 
related via q v (t — to) = —dS u (t — to)/dt. For the 
bimolecular reactions, q v {t — to) = <7 (t|r ,to) as 
given by Eg. l2*2l and S v (t — to) = S a (t\r ,to) as 
given by Eq. 1211 For the monomolecular reactions, 
?!/(*- t Q ) = qd(t\t ) = k d exp(-k d (t - t )) and 
S v (t - to) = exp(-fc d (t - t )). 

For each reaction R v , we generate a random num- 
ber £„, uniformly distributed in the interval < 
£i, < 1. If £i, < (1 — S u (oo)), a tentative next reac- 
tion time At„ is obtained from 

£„ = / q,(t')dt' = 1 - S v (At u ). (26) 
Jo 

If, however, £„ > (1 — S v (oo)), then the reaction 
R v does not occur and it is dropped from the list of 
possible reactions. From the remaining list of tenta- 
tive reactions, we choose as the actual next reaction 
the one that occurs first, 'provided that this reaction 
occurs within the maximum time step At max . Ac- 
cordingly, the system will be propagated through a 
time At as given by 

At = min({A^},Ai max ). (27) 

Note that if there is no tentative reaction for which 
the tentative next reaction time At„ < At max , then 
no reaction will occur within the time step. Here, 
we also mention for completeness that for associa- 
tion reactions S v (oo) ^ 0: for two particles, that 
can diffuse and react subject to the boundary con- 
dition Eq. 1191 there is a finite probability that they 
never react; this is related to the well-known fact 
that a random walker, that starts at the origin, can 
"escape" and never return to the origin. 

3. Propagate particles. The single particles are prop- 
agated according to pi(r, £|ro, to) in Eq.0 if a par- 
ticle decays, then the products are placed next to 
each other at positions centered around r. For each 
pair of particles, the following two substeps are ex- 
ecuted: 1) a new position for the coordinate R is 
obtained from Eq. El 2) if the pair has not re- 
acted, a new inter-particle vector r is obtained from 
p\{v, t|r , to) in Eq. |23 else, if it has reacted, the 
products are placed adjacent to each other at posi- 
tions around R. 

4. Update particles. Update identities of particles ac- 
cording to the executed reaction. Delete or add 
particles created or destroyed in that reaction. 



A proof of the algorithm is given in appendix B. 

Finally, we remark that the assumptions made above, 
namely that particles are spherical in shape and move 
by (free) diffusion, are not essential. An event-driven 
algorithm of this type could be set up in the case of 
non-spherical particles, interacting particles and/or those 
moving by other mechanisms than diffusion such as active 
transport. In these cases, the required Green's functions 
could be obtained numerically if necessary. 



III. RESULTS 

This section is organized as follows: first we study a 
simple bimolecular reaction to show that GFRD accu- 
rately reproduces theoretical results. Then we turn our 
attention to a very simple model of gene expression as 
a typical example of a system that is well handled by 
the GFRD technique. We specifically focus on the lev- 
els of noise in protein concentrations and find dramatic 
differences between GFRD and results from the chem- 
ical master equation, that ignores spatial fluctuations. 
Finally we compare the performance of GFRD to a con- 
ventional Brownian Dynamics algorithm. 



A. Bimolecular reaction 

To test the validity of our approach, we study the re- 
versible bimolecular reaction 

A + B^C (28) 

with forward rate constant k a and backward rate con- 
stant kd- As a first test of our algorithm, we focus on 
an isolated pair of particles A and B. We use a set up 
in which particle A, with diameter a, is placed at the 
origin and held fixed during the simulation. The second 
particle B, also with diameter a, is initially placed at 
random in a spherical shell of radius ro centered around 
particle A. We then propagate particle B for a time 
t s i m . During this time, particle B can diffuse freely with 
a diffusion constant D and it can associate with particle 
A with a rate constant k a and dissociate from it with 
a rate constant kd- Typically, particle B will associate 
with and dissociate from A a (large) number of times 
during the simulation. We repeat this whole procedure 
many times. This allows us to calculate the distribution 
function p Tev (r, t\ro, to), which gives the probability that 
the two particles A and B, separated by a distance ro at 
time to, are a distance r apart at a later time t. This 
numerical result can be compared to the analytical solu- 
tion recently derived by Kim and Shin for the reversible 
reaction shown in Eq. |2HI 10J . 

If the next reaction time were larger than the simula- 
tion time, t s im, then we could in principle directly propa- 
gate the particles through t s ; m . However, this would not 



7 



provide a stringent test of our algorithm. At each step, 
we therefore choose a maximum time step Ai max at ran- 
dom from the interval [10 _4 r, i s im]j where r = a 2 / D is 
the unit of time. This could be interpreted as mimicking 
the constraint on the maximum time step arising from 
the presence of other particles. 

Figure shows excellent agreement for p vev (r,t\ro,tQ) 
between theory and simulation for tq = a. We find sim- 
ilar agreement between theory and simulation for other 
initial distances tq and for other values of the diffusion 
constant D and reaction rates k a and k^. It should be 
realized that because the particles are initially placed 
at contact, many reactions can occur during the time 
t s i m . Moreover, because we divide the simulation time 
into smaller intervals, we must propagate the particles 
many times, using the Green's function for an extensive 
range of r, t — to and ro. Thus, at least for the case of 
an isolated pair of particles, this procedure provides a 
thorough test of our algorithm. 

Next, we want to study a more complex system in 
which a single particle of type A is held fixed at the cen- 
ter of a spherical container of radius R and is surrounded 
by Nb particles of type B. Particle A can again react 
with a particle B to form the product C according to the 
scheme in Eq. |2H1 particle C, like particle A, does not 
diffuse. Particles B and C do not react, although they 
are not allowed to overlap with each other. The excluded 
volume interactions between a pair of two B particles and 
between a pair of a B and a C particle is taken into ac- 
count by using reflecting boundary conditions, i.e. by 
setting k a — in Eq. |201 We note that the requirement 
that the B and C particles are not allowed to overlap, 



may impose a constraint on the maximum size of the 
time step, Ai max . The wall of the container is assumed 
to be reflecting. As no analytical solution exists for a 
pair of particles in the presence of a reflecting boundary, 
we introduce the further requirement that during a time 
step a particle can only interact with either the wall of 
the container or with another particle, but not with both. 

As the B particles diffuse through the container, they 
will come into contact with the fixed particle A. When in 
contact, the particles A and B can enter the bound state 
C with forward rate k a . When in the bound state, other 
B particles approaching the fixed C particle cannot react 
with it. Only after dissociation into the unbound state 
A + B, occurring at rate kd, can another reaction occur. 
On average, there will be a probability Abound of finding 
the A particle bound to the B particle. It is given by 

p b u ~ KNB (29) 

Abound- V + KNb ' 

where V is the volume available for the B particles, Nb 
is the total number of B particles and 

K = g AB {a)k a /k d (30) 

is the equilibrium constant. The function <?ab(0 is the 
radial distribution function for the pair of particles A and 
B. 

The radial distribution function gAB{f) describes the 
spatial correlations arising from the interactions between 
the particles K is conceivable that in this system 

the excluded volume interactions between the particles 
induce spatial correlations. These correlations could af- 
fect the density of B particles that are in contact with 



0.02 



0.01 




FIG. 2: The distribution p rcv (r, t\a, to = 0) for t = O.lr, It, 
10t and lOOr. For r < a, the distribution p Tev (r, t\cr, to) = 
due to the hard sphere repulsion between the particles. The 
bars denote the simulation results and the solid lines denote 
the analytical solutions of Kim and Shin [Tflj . Note that the 
particles are initially placed at contact (ro = cr). The forward 
rate constant k a = 1000 molecule - W -1 and the backward 
rate constant is kd = It . The unit of time r = a 2 /D. 
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FIG. 3: The probability Pbound that the A particle is bound to 
a B particle as a function of the total number of B particles for 
the reaction scheme shown in Eq. 1281 The symbols indicate 
the simulation results, while the solid line denotes the mean- 
field prediction (Eq. 1291 with <?a_b(<t) = 1). The radius of the 
container is R = 200<r and the equilibrium constant is chosen 
such that it is equal to the interaction volume V . The error 
bars in the simulation results are smaller than the size of the 
symbols. 



8 



the A particle and thereby the probability that the A 
particle is bound to a B particle. In Eq. [301 the distri- 
bution function at contact, gAB{&), thus describes the 
effect of the spatial correlations on the equilibrium con- 
stant. However, the concentrations that we consider here 
are very low and, as a result, the spatial correlations are 
expected to be small. Indeed, the simulations reveal that 
9ab{t) ~ 1 for all distances r. If <mb(c) = 1, then 
Eas. 1291 and 1301 reduce to the well-known mean-field re- 
sults that can straightforwardly be obtained by solving 
the macroscopic rate equations in steady-state. In Fig. [31 
we compare the simulation results to the mean-field pre- 
diction for Abound- We find excellent agreement. 

In conclusion, we have shown that our algorithm pro- 
vides an accurate way of simulating an assembly of par- 
ticles that can move by diffusion and react according to 
monomolecular and bimolecular reactions. As more com- 
plicated reactions, such as trimolecular reactions, can, in 
general, be decomposed into these elementary reactions, 
we are now in a position to simulate more complex sys- 
tems. 



B. Gene expression 

In this section we present results for a model of gene 
expression. It should be stressed that the model is highly 
simplified. The purpose here is to show the power of our 
approach. Nevertheless, we find interesting effects due to 
the spatial fluctuations of the components that could be 
of relevance for more realistic systems. 

The reaction network consists of the following reac- 
tions: 



A + B «=» C 

C -> P + A + B 

^prod 

P -> 

&dcc 



(31) 
(32) 
(33) 



In Eqs. 1311 1331 A represents a promoter region on the 
DNA and B a RNA polymerase molecule that moves by 
free diffusion and can bind with a forward rate k a to 
the promoter site to form the RNAp-DNA complex C. 
This complex can dissociate with a rate constant kd- Al- 
ternatively, it can produce a protein P at a production 
rate fc pro( j. Proteins degrade with a decay rate k df}c - Note 
that, in this model, when a protein is produced the RNAp 
dissociates from the DNA. 

In the living cell, the concentration of free RNAp - 
RNAp that is not bound to the DNA - is usually very 
low A s a result, spatial correlations are negligible 

(i.e., gAB{r) = 1) and the mean number of proteins, Np, 
can be obtained from the macroscopic rate equations. 
The result is: 



N P = K X K 2 



V + KxN B 



(34) 



Here K\ = k a /(k d + k pmd ) and K 2 = k prod /k dcc and N B 
is the total number of B molecules. 

In the simulations, we fix the promoter site, i.e. the A 
particle, in the center of a spherical container of radius 
R. The volume of the container is V — 1/im 3 , which is 
comparable to that of the Escherichia coli cell. The pro- 
moter site is surrounded by Np = 18 RNAp molecules, 
corresponding to the concentration of free RNAp of 30 
nM as found in the living cell ^| . The RNAp molecules 
move with a diffusion constant D = lfim s~ , which 
is comparable to that of similarly sized proteins |14| . 
We assume that, at contact, the RNAp can associate 
with the promoter site at a rate as determined by the 
Maxwell-Boltzmann velocity distribution. This leads to 
k a = it<j 2 (vab) = 3 • 10 9 M _1 s -1 , where vab is the rela- 
tive velocity of the particles A and B. We note that this 
estimate is equal to the rate of collisions between hard 
spheres in the low density limit 0, 0] . We could ar- 
rive at an alternative estimate for k a using the diffusion 
constant and a molecular "jump" distance A. This would 
lead to k a = 4tto- 2 D/\. Both estimates give similar re- 
sults for the value of k a . The dissociation rate is chosen 
such that the equilibrium constant K — k a /kd equals the 
one reported in |15|. yielding kd = 21.5s . We assume 
that the diameters of the promoter site and the RNAp 
molecules are equal and given by a = 5nm. 

Here, we only simulate the promoter site and the 
RNAp molecules explicitly in space. The proteins are 
assumed to be uniformly distributed in space. Moreover, 
we reduce both the degradation and the production of 
protein molecules to single-step Poisson processes. These 
assumptions are unrealistic. Nevertheless, this model al- 
lows us to demonstrate the power and the flexibility of 
our algorithm. In particular, the production and decay 
reactions can simply be added to our list of possible re- 
actions, (see section III Ef) . The next-reaction dis- 
tribution function for the production reaction is given by 
<7 P rod(i) = k prod N c exp(-kp ro dNct), where N c — if the 
RNAp is unbound and Nc = 1 if it is bound to the DNA, 
while the propensity function for the degradation reac- 
tion is given by q dcca , y (t) = k dccay N P exp(-k dccay N P t). 
In this way, the spatially-resolved GFRD scheme can nat- 
urally be combined with kinetic Monte Carlo schemes, 
such as the Gillespie algorithm |6|], that are based upon 
the spatially uniform chemical master equation. 

In figure 0] we show the mean number of proteins Np 
as a function of the protein production rate fc pro d, while 
keeping the decay rate fixed at fcdecay = 0.04s -1 . As the 
concentration of the RNAp is low and spatial correlations 
are expected to be negligible, the simulation results for 
the average number of proteins, Np, should follow the 
mean- field prediction of Eq.021 Fig. 0] shows that this is 
indeed the case. However, in contrast to the mean-field 
analysis, the GFRD simulations allow us to quantify the 
effect of the spatial fluctuations of the RNAp molecules 
on the noise in the protein synthesis. 

We can quantify the magnitude of the noise in protein 
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production by computing the following quantity |15| : 



N 2 P {t) - Np 
Np 2 



(35) 



In the following analysis we have changed the degradation 
rate such that the average number of proteins, Np, is 
constant. This allows us to focus on the effect of spatial 
fluctuations on the noise in protein production - we thus 
eliminate the fairly trivial changes in the noise due to 
changes in the average number of proteins. 

Since we are interested in the importance of spatial 
fluctuations in gene expression, it is natural to compare 
the GFRD results to those obtained using the chem- 
ical master equation. The latter approach does take 
into account the discrete nature of the reactants and the 
stochastic character of their interactions, but it treats 
the spatial fluctuations in a mean-field manner: at each 
instant, it is implicitly assumed that the particles are uni- 
formly distributed in space. This approach is justified if 
there are many non-reactive collisions to stir the system 
in between the reactive collisions. However, the RNAp 
is present in low copy numbers, and, upon contact, it 
rapidly associates with the promoter site on the DNA. 
As a consequence, this reaction is diffusion-limited. This 
could have important implications for the noise in gene 
expression. 

Using the techniques described in Q, we can analyti- 
cally obtain the noise from the chemical master equation. 
It is given by 



VP 



Cprod 



k„.N K 



Np k pmd k a N B + N P (k a N B + k d + k prod ) 



(36) 

The first term on the right describes the result that would 
have been obtained if gene expression were a simple lin- 
ear birth-and-death process. The second term reflects the 
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FIG. 4: The mean protein number Np as a function of the pro- 
tein production rate fc pro d as obtained from the GFRD sim- 
ulations for the reaction scheme shown in Eqs. 1311 1331 The 
solid line denotes the mean-field prediction given by Eq. 



fact that in order to produce a protein, it is necessary, 
albeit not sufficient, for a RNAp molecule to bind to the 
promoter site. This term, and thus the noise in gene ex- 
pression, goes through a minimum at k pTO d = k a Np + kd 
and vanishes for both small and large k pro d- I n these 
regimes, gene expression reduces to a simple linear birth- 
and-death process. In the limit of small fc pro( j, the pro- 
duction of the protein is the rate limiting step. The 
RNAp molecule will associate to and dissociate from the 
promoter site a large number of times, before it actually 
induces gene expression. The former process is thus in 
equilibrium on the time scale of gene expression. Hence, 
the birth term is given by /c b irth = Pbound&prod, with 
Pbound being the probability that a RNAp molecule is 
bound to the promoter (see Eq. I29[) : the death term is 
given by fcdeath = fcdocay In the limit of large k pmA , the 
binding of a RNAp molecule to the promoter site is the 
rate limiting step: as soon as a RNAp molecule is bound 
to the promoter, a protein will be produced. This means 
that the birth term is given by fcbirth = k a (l — Pbound)[B]; 
the death term is again fcdeath = fcdecay For a linear 
birth-and-death process, the noise is determined by the 
average number of proteins, r\p = \j \J Np 3J. As we 
have set the decay rate fcdecay such that Np is constant, 
the noise in gene expression must be the same in the 
limiting regimes of small and large fc pro d, in which gene 
expression reduces to a birth-and-death process. 

Figure [5] shows the noise in the protein concentration 
as a function of the synthesis rate for Np = 1000. The 
GFRD results are compared to those obtained using the 
chemical master equation. It is seen that for small fc pro d 
both approaches yield identical results. In this regime, 



GFRD, D=0.1 urrT % 



GFRD, D=1 urrT s~ 



Master equation 



FIG. 5: The noise in protein level r\p as a function of protein 
production rate fc pro d for the reaction scheme shown in Eqs. 1311 
- 1331 Compared are the results obtained by GFRD with a 
diffusion constant of D — l/im 2 s _1 (■) and D = 0.1/im 2 s -1 
(*) and the result using the chemical master equation (dashed 
line). The mean number of proteins was held constant at 
Np — 1000 by changing the degradation rate fcdecay of the 
protein. 
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protein synthesis is the rate-limiting step. Indeed, on the 
time scale of gene expression the RNAp molecules have 
sufficient time to become well mixed in the cell. As a 
result, the effects of diffusion are negligible and the noise 
reduces to the expected value for a linear birth-and-death 
process. 

However, for k pro d > Is -1 , spatial fluctuations can 
have a drastic effect on the noise in gene expression. In 
this regime, the noise of the spatially-resolved model is 
larger than that of the "well-stirred reactor" model. In 
fact, Fig. [S] shows that it grows fairly rapidly with in- 
creasing fcp r od- The increase in noise is due to a very 
broad distribution of arrival times of RNAp molecules at 
the promoter site, much broader than the corresponding 
Poisson distribution for the system without spatial fluc- 
tuations. It is also seen that for very large production 
rates, the noise ultimately reaches a plateau value. At 
these high values of fc pro d the promoter site becomes an 
"absorbing" boundary for the RNAp molecules. Fig. 
also reveals that the height of the plateau increases as the 
diffusion constant D becomes smaller. This is not sur- 
prising, because the importance of spatial fluctuations 
is expected to be larger for smaller diffusion constants. 
However, it does clearly show that in order to determine 
the significance of spatial fluctuations in gene expression, 
it is of interest to establish the value of the diffusion con- 
stant of the RNA polymerase experimentally. 

This model of gene expression is obviously highly sim- 
plified - our aim here is to present a new particle-based 
technique to simulate biochemical networks. Neverthe- 
less, the results show that spatial fluctuations are poten- 
tially important in gene expression. Moreover, GFRD 
will now make it possible to study systematically how 
significant the effects of spatial fluctuations are on the 
noise in gene expression using more refined models. 



C. Performance 

The essence of the GFRD scheme is to exploit the ana- 
lytical solution of the Smoluchoswki equation for a pair of 
interacting particles to set up an event-driven algorithm. 
This allows GFRD to make large jumps in time when 
the particles are far apart from each other. Clearly, the 
performance of the algorithm depends on the density of 
the system: the further the particles are apart, the larger 
the time step that can be used and the better GFRD will 
perform in comparison to brute-force Brownian Dynam- 
ics schemes (see section III All . It is thus of interest to 
compare the distribution of propagation times in GFRD 
to the time step used in a brute-force Brownian Dynam- 
ics scheme as a function of density. 

In figure we show the distribution of propagation 
times At for the bimolecular reaction described in section 
iHlAl as a function of density. For N B = 18 ([£?] = 30 
nM), the value used in the above models of gene expres- 
sion, the distribution has a maximum at At = 1 • 10 _4 s. 
The propagation times become smaller if the density in- 
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FIG. 6: The distribution of propagation times At for a system 
consisting of a single particle A in the center of a spherical 
box of volume V = 1/xm 3 , surrounded by Nb particles B; 
the particles A and B can react according to a bimolecular 
reaction scheme (see Eg. 1281 . The distributions are shown for 
N B = 18 (bold solid line), N B = 50 (solid line) N B = 125 
(dashed line), N B = 300 (dashed-dotted line) and N B = 600 
(dotted line), corresponding to a concentration of [B] = 30 
nM, [B] = 83 nM, [B] = 210 nM, [B] = 500 nM, and [B] = 
1000 nM, respectively. In the GFRD simulations, we use a 
lower cut-off for the propagation time, 2.5 ■ 10~ 10 s, which 
corresponds to the time step used in a brute-force Brownian 
Dynamics simulation. 



creases. For [B] = 1/xM, the peak in the distribution 
shifts to At = 1 • 10" 6 s. 

In biochemical networks, the concentrations of the 
components can be very low. In gene networks, for ex- 
ample, the concentrations of the gene regulatory proteins 
are often in the nM regime. In signal transduction net- 
works, the concentrations of the components may also be 
fairly low, i.e. in the /xM range. The analysis presented 
here suggests that with GFRD it should be possible to 
reach time steps of at least 10~ 6 — 10 _4 s for such net- 
works. In contrast, in a brute-force Brownian Dynam- 
ics simulation, we cannot use a time step larger than 
lO" 10 - 10~ 9 s (10~ 5 - lQ- 4 a 2 /D) in order to preserve 
the correct distribution (as, for instance, defined by the 
requirement that the analytical solution for the Green's 
function p rev (r, t\ro, to) as shown in Fig. [SJcan be accu- 
rately reproduced). We thus believe that under biologi- 
cally relevant conditions, GFRD can be three to six or- 
ders of magnitude faster than conventional particle-based 
schemes for simulating biochemical networks in time and 
space. 



IV. CONCLUSIONS 

We have developed a new technique, called Green's 
Function Reaction Dynamics, to simulate biochemical 
networks at the particle level and in both time and space. 
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The main idea of the technique is to choose a maximum 
time step such that only single particles or pairs of parti- 
cles have to be considered. For these particles, the Smolu- 
chowski equation can be solved analytically using Green's 
functions. The analytical solution can then be used to 
set up an event-driven algorithm, quite analogous to the 
kinetic Monte Carlo schemes as originally developed by 
Bortz, Kalos and Lebowitz | 19| to simulate Ising spin 
systems and by Gillespie to numerically solve the chemi- 
cal master equation We would like to stress, however, 
that in contrast to the widely used "Gillespie" algorithm, 
our technique makes it possible to simulate biochemical 
networks at the particle level and in both time and space. 

The analysis presented in section IIII CI shows that 
GFRD is highly efficient. This should make it possible 
to simulate biochemical networks at much larger length 
and time scales than hitherto possible. In addition, we 
believe that the scheme has the potential to be even more 
efficient. In the current scheme, we use a global maxi- 
mum time step that pertains to all particles in the sys- 
tem. It seems natural, however, to assign to each parti- 
cle an individual maximum time step. In such a scheme, 
each particle would have its own individual clock. This 
approach would make it possible to devote most compu- 
tational effort to those particles that interact frequently; 
the particles that are initially far from other particles are 
only updated once the time has come when they have 
a reasonable chance to interact. A second possible im- 
provement would be to exploit the low concentration of 
the components in another way. In the current scheme, 
we explicitly take into account excluded volume interac- 
tions. In fact, this often imposes a limit on the maxi- 

I 



and where P n is the nth Legendre polynomial, J n and 
Y n are the nth Bessel function of the first and the second 
kind, D = Da + F>b is the total diffusion constant of the 
two particles A and B and a = (cU + de)/2, where d,A 
and ds are the diameters of the particles A and B, re- 
spectively; here, and below, we take to = 0. The Green's 
function can be expressed in a more compact notation by 



oo 

P r 2 (r,d,ct>,t\r ) = Y / CnPn{cos6)R n (r,t). (39) 

n=0 



mum time step. If the concentrations are low, however, 
we would expect the excluded volume effects to be neg- 
ligible for the behavior of the network. In future work, 
we will show how these observations can be incorporated 
into the algorithm to even further enhance the perfor- 
mance of Green's Function Reaction Dynamics. 
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Appendix A: Solution of Smoluchowski equation 

Here we consider a pair of particles A and B that move 
by free diffusion, but, upon contact, can react with a rate 
constant k a . The Green's function p 2 (r, t|ro, to) for this 
pair can be obtained by exploiting the analogy between 
the diffusion of particles and the conduction of heat in 
solids. The corresponding Green's function for the con- 
duction of heat in solids is derived in [llj . The Green's 
function p 2 (r, t|ro, to) is 



(37) 



(38) 



The probability f(r\r ,t) of finding the particles sepa- 
rated by a distance between r and r + dr at time t is 
given by 

oo 

f(r\r ,t) = 2vr^a i g„(7r)r 2 i? n (r,t), (40) 

n=Q 

with 

Qn{0)=( sin 6P n (cos 6)d6. (41) 
Jo 

The conditional probability g{9\r, r$, t) that two particles 
are at an angle between 9 and 9 + d9 with respect to 



p r 2 {r,9,0,t\r o ) = VY2n + l)P w (cosfl) / e - Du2t F n+1/2 (ur)F n+1/2 (ur )udu, 

^V^t^ Jo 

where 

{2ak a + l)[J v (ur)Y v (w) - Y v {ur)J u {ucr)] - 2ua[J u {ur)Yl{ua) - Y v {ur)J' v (w)] 



F v {ur) 



([(2ak a + 1) J u (ua) - 2uaJ' u {ua)] 2 + [{2ak a + l)Y v {ua) - 2uaYl(ucr)} 2 ) 1 / 2 

I 
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the original direction i"o = yb — ya, given that they are 
separated by a distance r at time t, is 

oo 

g(6\r,r ,t) = 2tt £ C n Q' n (9)r 2 R n (r, t). (42) 

71=0 

The survival probability S(t) is given by 

/>oo 

= / f(r\r ,t)dr. (43) 

The above integral is complicated but it follows from the 
properties of the diffusion equation that it must be identi- 
cal to the familiar survival probability for the spherically 
symmetric case 0, . 

The above distribution functions suggest a straightfor- 
ward procedure for drawing a new position r from the 
Green's function p2( r j*l r o)- We pretabulate Q n {9) and 
R n (r, ro , t) up to a certain order N. From this we con- 
struct the probability distribution f(r\ro,t) and we draw 

I 



where 

Rl = (2crk a + l)J n+ i(ua)-2u(jJ' T ^i(ua), (47) 
R 2 = (2ak a + l)Y nh i(ua) - 2uaY^_(ua), (48) 

^1 = J n+i ( ur )Jn+± ( ur o) - Y n+± ( ur ) Y n+± (Wo), (49) 

F 2 = J„ + i (wr)y^ (ur ) + J n+ i (ur Q )Y n+ i (ur) . (50) 

As the correction term p cor r is usually rather small for 
small t, the number of terms N that needs to be in- 
cluded in order to compute the functions f(r, \ro,t) and 
g(9\r, ro, i), is strongly reduced. 

Appendix B: Proof of the algorithm 

Let P(t, n)dt be the probability that the next reaction 
will occur in the time interval between t and t + dt and 
will be reaction R^. As described in section lTl Bl the time 



a new distance r from this distribution. Next, we draw 
9 from the distribution g(6\r,ro, t) and finally we choose 
(j) uniformly distributed between and 2tt. 

This procedure works well for large t. For small t, how- 
ever, the above procedure becomes rather cumbersome as 
the number of terms N that needs to be included in order 
for the summations in Eas. 1401 and 1421 to converge, be- 
comes very large. The reason is that ( r i *l r o) becomes 
a sharply peaked function around ro for small t. How- 
ever, for small t, the probability that the two particles 
will interact with each other, is relatively small. In other 
words, for small t the full solution ^(r, t|ro), is domi- 
nated by free diffusion. We can exploit this observation 
in order to reduce N by writing the Greens's function as 
(r,t|ro)+Pcorr(r,£|ro), where pf roc is the 
solution for free diffusion and p C orr is a correction term 
that takes into account the reacting boundary at r = a. 
The free diffusion term can easily be computed from 



(44) 



(45) 



(46) 



step is chosen such that the reactions occur independently 
from each other. The probability Pit, [i)dt is therefore 
given by 

M 

P(t, l x)dt = q ll {t)dtJ[S v (t). (51) 

v=l 

The algorithm features a maximum time step. It is thus 
conceivable that not a single reaction occurs within a 
time step. The probability Q(Ai max ) that no reaction 
occurs within a time step of size At max is given by 

M 

Q(At max ) = J] S„(Ai maJ£ ). (52) 

!/=l 

It can easily be shown that the procedure outlined 
in section III El is consistent with Eqs. and [521 Let 
<2(Ai max ) be the probability that the procedure de- 
scribed above does not yield a reaction within the time 



Pbee(r, 9, 4>, t\ro) 



■ exp 



r2 + r o ~ 2rr cos ( 



4Dt 



[AnDtf /2 

Using the fact, that pf ree can also be written as 

Ph ee (r, 9, 4>, t\r ) = -=y^(2n+ l)P„(cos0) / e~ Du * J n+ i/ 2 (w) J n+ i/ 2 (ur )udu, 

we find by comparison with Eq. [37| that p cor r can be expressed as 



Peon (r, 9, 4>,t\r q) = - 



4tt 



-— y^{2n + l)P n (cos9) / 

A^orrr, Jo 



Rl + Rl 



[RlFx + R 2 F 2 )udu, 
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step of size Ai max . It is given by 

Prob (t v > At max for all v) (53) 

M 

[] Prob > (1 - S^AU))) (54) 

M 

Y[ S v (At max ) (55) 

Q(Ai max ). (56) 

Similarly, let P(t, [i)dt be the probability that the above 
described procedure yields, at time t, reaction as the 
next reaction . It is given by 

P(t, n)dt = Prob (t < < t + dt) 

x Prob (t„ > £ M for all i/ ^ /i) (57) 

M 

= <Z„(t)dt II Prob & > ( X - ^(^))X58) 

iy = l 
M 

= g/i (t)dt n ^(^) (59) 

v = \ 

= P(t,n)dt. (60) 



Q(At max ) = 
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